#! /usr/bin/env python

# NMDB stations: get current coordinates, cut-off, ...
# Copyright (C) 2010-2011 Christian T. Steigies <steigies@physik.uni-kiel.de>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

## --------------------------------------------------------------------------

#http://mysql-python.sourceforge.net/MySQLdb.html
import MySQLdb # python-mysqldb
import Gnuplot # python-gnuplot
import os
import optparse
import ConfigParser

from datetime import datetime, timedelta
from numpy import * # python-numpy

## --------------------------------------------------------------------------

def config_section_map(section,Config):

    dict1 = {}
    options = Config.options(section)
    for option in options:
        try:
            dict1[option] = Config.get(section, option)
            if dict1[option] == -1:
                DebugPrint("skip: %s" % option)
        except:
            print("exception on %s!" % option)
            dict1[option] = None
    return dict1

## --------------------------------------------------------------------------

def parse_options():

    parser = optparse.OptionParser("usage: %prog [options] [dummy]")
    parser.add_option("-H", "--host", dest="hostname", default="db02.nmdb.eu", type="string", help="specify hostname to run on")
    parser.add_option("-P", "--port", dest="portnum", default=3306, type="int", help="port number to run on")
    parser.add_option("-f", "--file", dest="rcfile", default=os.environ["HOME"]+"/.nmdbrc", type="string", help="config file to use")
    parser.add_option("-d", "--date", dest="date", default="", type="string", help="show stations with start_date before this date")
    parser.add_option("-l", "--label", action="store_true", dest="label", default=False, help="also plot station labels")
    parser.add_option("-c", "--coast", dest="coast", default=1, type="int", help="coast type")
    parser.add_option("-p", "--projection", dest="projection", default=1, type="int", help="projection type")
    parser.add_option("-i", "--isolines", action="store_true", dest="isolines", default=False, help="plot cutoff isolines")
    parser.add_option("-C", "--contours", action="store_true", dest="contours", default=False, help="plot cutoff contours")
    parser.add_option("-v", "--verbose", action="store_true", dest="verbose", default=False, help="make lots of noise")
    
    (options, args) = parser.parse_args()

    return (options)

## --------------------------------------------------------------------------

def parse_config(rcfile):

    Config = ConfigParser.ConfigParser()
    Config.read(rcfile)

    username = config_section_map("alert",Config)['username']
    password = config_section_map("alert",Config)['password']

    return(username,password)

## --------------------------------------------------------------------------

(options) = parse_options()
(username,password) = parse_config(options.rcfile)


db=MySQLdb.connect(host=options.hostname,port=options.portnum,db="nmdb",user=username,passwd=password)
c=db.cursor()

## --------------------------------------------------------------------------
##describe station_information;
##+-----------------------------+---------------+------+-----+---------+----------------+
##| Field                       | Type          | Null | Key | Default | Extra          |
##+-----------------------------+---------------+------+-----+---------+----------------+
##| ID                          | bigint(20)    | NO   | PRI | NULL    | auto_increment |
##| station_full_name           | varchar(100)  | NO   |     |         |                |
##| station_short_name          | varchar(20)   | NO   |     |         |                |
##| start_date_time             | date          | NO   |     | NULL    |                |
##| end_date_time               | date          | YES  |     | NULL    |                |
##| head_organization           | varchar(1000) | NO   |     |         |                |
##| principal_investigator      | varchar(200)  | NO   |     |         |                |
##| contact_person              | varchar(200)  | NO   |     |         |                |
##| latitude_deg                | float         | NO   |     | 0       |                |
##| longitude_deg               | float         | NO   |     | 0       |                |
##| altitude_m_asl              | float         | NO   |     | 0       |                |
##| geomag_cutoff_GV            | float         | NO   |     | 0       |                |
##| avg_count_per_sec_solar_min | float         | NO   |     | 0       |                |
##| avg_count_per_sec_solar_max | float         | NO   |     | 0       |                |
##| ref_pressure_mbar           | float         | NO   |     | 0       |                |
##| barometric_coefficient      | float         | NO   |     | 0       |                |
##| detector_housing            | varchar(1000) | YES  |     |         |                |
##| counter_tube_type           | varchar(100)  | NO   |     |         |                |
##| counter_tube_num            | int(11)       | NO   |     | 1       |                |
##| detector_info               | varchar(1000) | NO   |     |         |                |
##| pressure_sensor             | varchar(1000) | NO   |     |         |                |
##| clock_info                  | varchar(200)  | NO   |     |         |                |
##| additional_station_info     | varchar(2000) | NO   |     |         |                |
##+-----------------------------+---------------+------+-----+---------+----------------+
##23 rows in set (0.00 sec)

if options.date:
    today = options.date
else:
    today = datetime.utcnow()
print "NMDB stations at", today, "\n"

# get list of stations with station information
querystring = "SELECT station_short_name FROM station_information GROUP BY station_short_name;"
c.execute(querystring)
station=[]
answer=True
while answer:
  answer=c.fetchone()
  if answer:
    station.append(answer[0])

## --------------------------------------------------------------------------
dat = open('nmdb_stations.d','w')
loc = open('nmdb_stations.loc','w')

## get current coordinates for each station
for i in station:
    QUERYSTRING = ("SELECT start_date_time, latitude_deg, longitude_deg, "
                   "altitude_m_asl, geomag_cutoff_GV, station_short_name "
                   "FROM station_information "
                   "WHERE start_date_time <= \"%s\" "
                   "AND station_short_name ='%s' "
                   "ORDER BY start_date_time DESC LIMIT 1;") % (today, i)
    if options.verbose:
        print QUERYSTRING
    c.execute(QUERYSTRING)
    answer=c.fetchone()
    if answer:
        start_date = answer[0]
        latitude = answer[1]
        longitude = answer[2]
        altitude = answer[3]
        cutoff = answer[4]
        name = answer[5]
        if options.verbose:
            print i, ":", name, latitude, longitude, altitude, cutoff

        ## x, y, z
        print >>dat, longitude, latitude, 0

        ## x, y, size, angle, fontno, justify, text
        print >>loc, longitude, latitude, 5, 0, 5, 0, name

## --------------------------------------------------------------------------
## create GMT script
out = open('nmdb_stations.sh','w')
print >>out, """#!/bin/sh
# NMDB station map

PATH="/usr/lib/gmt/bin:${PATH}"
AWK=${AWK:-awk}

FILE="nmdb_stations"

# map area
RANGE="-R-180/180/-90/90"
"""
if options.coast == 1:
    print >>out, 'COAST="-Dc -G160 -A10000"'
elif options.coast == 2:
    print >>out, 'COAST="-Dc -B45g45/30g30 -A10000"'
else:
    print >>out, 'COAST="-Dc -G160 -A10000"'

if options.projection == 1:
    print >>out, 'PROJECTION="-JX6id/3.75id"'
elif options.projection == 2:
    print >>out, 'PROJECTION="-Jj0/0.025i" # Miller'
elif options.projection == 3:
    print >>out, 'PROJECTION="-JR0/8i" # Winkel'
elif options.projection == 4:
    print >>out, 'PROJECTION="-JN0/8i" # Robinson'
else:
    print >>out, 'PROJECTION="-Jx6id/3.75id"'

print >>out, """VIEW="-E180/90"

gmtset PAGE_ORIENTATION portrait
gmtset PAPER_MEDIA A4+
gmtset PLOT_DEGREE_FORMAT ddd:mm:ss
gmtset BASEMAP_AXES WESN
gmtset D_FORMAT %3.0f
"""
if options.contours or options.isolines:
    print >>out, """# calculated with MAGNETOCOSMICS
${AWK} '{if (NR>5) print $2*1, $1*1, $4*1}' CutoffVsPosition.out > cutoff.d

# cutoff
z0=`minmax cutoff.d -C -I25 | $AWK '{print $5}'`
z1=`minmax cutoff.d -C -I25 | $AWK '{print $6}'`

makecpt -Cno_green -T$z0/$z1/1 > topo.cpt
"""
if options.contours or options.isolines:

    print >>out, """ # paint surface (-I)
#pscontour ${RANGE} ${PROJECTION} cutoff.d -Ctopo.cpt -I -K  > ${FILE}.ps
# overlay contours (no -I, set linewidth with -W)
#pscontour ${RANGE} ${PROJECTION} cutoff.d -W0.1p -Ctopo.cpt -K -O >> ${FILE}.ps
pscontour ${RANGE} ${PROJECTION} cutoff.d -W0.1p -Ctopo.cpt -K > ${FILE}.ps

# continent outlines
pscoast ${RANGE} ${VIEW} ${PROJECTION} ${COAST} -W0.25p -P -K -O >> ${FILE}.ps

# overlay station locations
psxyz -G102/0/255 nmdb_stations.d -W0.5p -Ss0.15 ${PROJECTION} -R -K -O >> ${FILE}.ps

# overlay station names
pstext -S0.5p -G255/0/0 nmdb_stations.loc -R ${PROJECTION} -O >> ${FILE}.ps

"""

## --------------------------------------------------------------------------
  
c.close()

## --------------------------------------------------------------------------
